Snow metamorphism: a fractal approach 
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Abstract 

Snow is a porous disordered medium consisting of air and three water phases: ice, vapour 
and liquid. The ice phase consists of an assemblage of grains, ice matrix, initially arranged over 
a random load bearing skeleton. The quantitative relationship between density and morpho- 
logical characteristics of different snow microstructures is still an open issue. In this work, a 
three-dimensional fractal description of density corresponding to different snow microstructure 
is put forward. First, snow density is simulated in terms of a generalized Menger sponge 
model. Then, a fully three-dimensional compact stochastic fractal model is adopted. The 
latter approach yields a quantitative map of the randomness of the snow texture, which is 
described as a three-dimensional fractional Brownian field with the Hurst exponent H varying 
as continuous parameters. The Hurst exponent is found to be strongly dependent on snow 
morphology and density. The approach might be applied to all those cases where the morpho- 
logical evolution of snow cover or ice sheets should be conveniently described at a quantitative level. 
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I. INTRODUCTION 



Physical and mechanical properties of snow greatly vary over space and rapidly evolve in 
time (snow metamorphism), affecting thermal conductivity, strength, heat capacity, density 
and texture [IH3]. Density and texture depend on meteorological conditions, temperature, 
pressure and humidity, which, together with the mechanical strains cause metamorphism 
and, ultimately, instability of snow [4-14J. Density is the parameter used for classifying 
and assessing snow properties, thanks to the simplicity of in-situ measurements and can be 
adopted as a parameter for quantifying characteristics as viscosity, shear stress and strength, 
cohesion and mechanical properties, such as Young's Modulus and Poisson's ratio. Usually, 
researchers refer to the specific density, defined as the ratio of the snow density p sn0 w to 
the ice density pj ce =917 kg/m 3 and ranging between 0.05 to 0.60, to describe mechanical 
properties of snow. However, different types of snow having the same density might exhibit 
completely different mechanical features. Elucidating the quantitative relation between den- 
sity and mechanical characteristic of snow is an open issue worthy of investigation beyond 
purely speculative interests. For example, the knowledge of density of added or lost snow 
is needed in altimeter measurements of dynamic thinning of Arctic and Antarctic ice sheets 
and monitoring of incipient avalanches [TBHTTj . 

Scaling properties and fractional calculus have been extensively adopted to characterize 
different classes of materials [IB]. In particular, fractal concepts have been already used 
for modelling snow crystals and reproducing snowflake morphology (see e.g. the von Koch 
snowflake curves and the application of general iterated function systems). In the last two 
decades, fractal theories of snow have resulted in several applications: from the measure 
of fractal dimension of images, to remote sensing and mapping of snow cover and depth 
distribution by satellite and Lidar images, from the study of the roughness of the snow- 
pack, to the determination of the air flux across snow surface and to the definition of the 
spatial variation of the snow water equivalent [THH2E]. The fractal character of snow has 
been recently studied by the Centre Etudes de Neige of Meteo France [271 EE] by analysing 
three-dimensional tomography obtained for cubic samples with 2.5 mm size and different 
densities. Relationships between mechanical properties (tension, shear strength, toughness) 
and specific density in terms of power laws with non integer exponents suggesting fractal 
features of snow are reported in J2HH3I]- Furthermore, fractal theory has been used to inves- 
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tigate the stability of the snow cover at larger scales. Based on the renormalisation group 
model, the probability of occurrence of snow avalanche events, assuming scale invariance of 
the snowpack at the smallest scales with a consequent implication of the same behavior at 
the largest scales [32j 133] . 

The present work is addressed to investigate the multiscale character of snow density, by 
adopting a fractal description of the distribution of ice grains, able to reproduce the local 
randomness of real microstructure. The paper is organized as follows. In Section II, snow 
density is simulated by means of a generalized Menger sponge model, characterized by a 
discrete set of Hurst exponent values. In Section III, a fully three-dimensional fractional 
Brownian model is reported. This model has the advantage to reproduce the randomness of 
the local microstructure of snow samples in a more realistic way. Snow density is mapped 
to the three-dimensional fractional Brownian field, with the Hurst exponent H continuously 
ranging from to 1. The proposed model provides a fully three-dimensional analysis in terms 
of a continuum fractional Brownian field rather than the discrete Menger model description 
[SUES]- I n the framework of this model, fractal dimension D and Hurst exponent H quantify 
ice distribution and reproduce the values of snow density for different microstructures. This 
fully three-dimensional fractal model can be used to relate snow texture, obtained by in-situ 
measurements o remote imaging techniques, and density changes, important for monitoring 
snow cover and ice sheets dynamics, through the estimate of Hurst exponent at different 
snow metamorphism stages. 



II. MENGER SPONGE MODEL 



In this Section, the fractal, known as Menger sponge, is used to characterize scale invariant 
features of porosity and density of snow. The Menger sponge is generated as shown in the 
scheme of Fig.l. First, the bulky cube (a) is divided into 3 x 3 x 3 = 27 equal subcubes (b). 
Then, 7 subcubes are removed from the center of each face and from the center of the cube, 
resulting in Nf = 20 filled subcubes and N e = 3 3 — Nf = 7 empty subcubes. This single 
step is repeatedly applied to the remaining cubes. For a solid cube with linear size r , the 
first step Menger sponge is characterized by linear size T\ = 3r (Fig. [ljb)) . The second step 
Menger sponge is characterized by linear size r2 = 9tq as shown in Fig. [l](c). In general, 
i~i = (3*)ro is the linear dimension of the fractal cube (Menger sponge) at the iteration i. 
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The fractal dimension of the Menger sponge is given by D = In Nf/ In 3 ~ 2.727. Here, 
we assume self-similarity which implies a linear relationship H = 3 — D ~ 0.273 between 
fractal dimension D and Hurst coefficient H. 

The Menger sponge has been widely used to model porous media, whose relevant param- 
eter is the porosity [18] . The porosity <fii, defined as the relative volume of voids per unit 
volume, can be expressed by the following relationship: 



(1) 



Pi 

where po is the initial density of the bulky cube and pi is the density of the Menger sponge 
at the iteration %. By taking into account that pi and po are inversely proportional to the 
volumes and, then, to the linear sizes tq and r*, Eq. can be written as: 

* = !-(*) , (2) 

where D is the fractal dimension, that is the scale-independent parameter characterizing the 
morphology of a porous material. Therefore, density and void index of the Menger sponge, 
as a function of the linear dimensions of the solid cube ro and fractal cube r\ can be written 
at the i th iteration as: 

In N f 

For the Menger sponge shown in Fig. pjb), the porosity is <fri = 7/27 and the density is 
pi = 20po/27, while for the Menger sponge of Fig. [jjc), the porosity is 02 = 329/729 and 
the density is p 2 = 400p /729. 

The above described procedure can be generalized by removing an arbitrary number N e 
of subcubes (instead of 7) out of an arbitrary number of solid cubes Nf (instead of 20). This 
generalized construction results in fractal structures with Hurst exponent different than 
H = 0.273. 

To model snow samples, an homogeneous ice cube characterized by density po = Pice = 
917 kg/m 3 and linear size r = r grain is considered at the initial step. Then, snow is obtained 
as a Menger sponge, i.e. a fractal form of ice, characterized by density p S now = Pit 

and linear 

size r snow = Vi at the iteration i. Thus, Eq. (J3]) is rewritten for snow as follows: 



In Nf 

, 3--^f / \ 3-D / \ H 

ysnow I 1 gram \ j 1 gram \ j 1 gram 

Pgrain 



(4) 



Fractal dimension D and Hurst exponent H can be taken as a measure of porosity of the 
snow sample. The fractal character of snow can be analysed by using images obtained on 
various cubic samples with same size and different densities, showing the granular structure 
and the spatial distribution of its voids. By means of the box-counting method, the fractal 
dimension of four snow samples, characterized by grains with different diameter, could be 
determined [28J. By applying the generalized Menger sponge model the fractal dimension 
D of different classes of snow can be calculated (Table I). These values have been obtained 
by using the same size of the samples: r snow = n = 100 mm and the grain size r grain ranging 
between 0.05 mm and 0.25 mm. By analyzing density and void index as a function of the 
linear sample dimension rj (grain size), we observe that as the grain size increases, snow 
differs more and more from ice. We also observe that the values of the fractal dimension 
D measured by the box-counting method, ranging between D = 2.62 to D = 3 [28J, are 
consistent with the values calculated by the Menger sponge model (Table I). At small scales, 
ice and snow approximately show the same behavior while the spatial variability of the 
density does not greatly influence the mechanical properties. Therefore, we argue that snow 
density is a function of the scale and the probability to find large defects (e.g. super- 
weak zones in a weak layer) increases with the dimension of the snow grain, as for example 
providing more intrinsic brittleness for large snow slopes J32j[33]. Numerical results reported 
in Table I confirm that D is an accurate measure of the distribution of the ice mass into 
snow samples. 



III. THREE-DIMENSIONAL FRACTIONAL BROWNIAN MODEL 

In the previous section, the generalized Menger sponge model, characterized by scale 
invariant porosity, is used to develop a fractal description of snow density of samples with 
different microstructures. This description is of significant practical relevance, however it 
does not fully capture the fractal structure of snow samples. One drawback is due to the 
discreteness of fractal dimension. Moreover, the porosity of a real-world fractal should be 
free of scale requiring an infinite number of iterations for generating the sponge. On account 
of these limitations, in order to describe snow as a sintered porous material consisting of 
a continuous ice network, a generalization of the random midpoint displacement algorithm 
is implemented to model snow as a three-dimensional fractal heterogeneous medium [35J. 
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The method for generating compact fractal disordered media relies on fractional Brownian 
functions, which are characterized by correlation depending on the distance r as a power 
law. The approach is based on the function: 

f H ,j ( r ) = ^dYl ^ + a i' d > ( 5 ) 

k 

with r = (iijiijis). The sum is calculated over the k endpoints of the lattice and the 
quantity is a random variable defined at each iteration j as: 

°l d = °l{^) [1"2 2( ^], (6) 

where the quantity o"o is drawn from a Gaussian distribution with zero mean and unitary 
variance. The Hurst exponent H, ranging from to 1, is the input of the algorithm which is 
implemented at each iteration j according to the following procedure. Initially, the lattice 
is fully homogeneous, with the function describing the fractal property taken as a 

constant, e.g. fu(j') = 0. Then, the values of the function fujij') are seeded as random 
variables at the eight vertices of the cube. The value assigned to the central point is obtained 
by means of Eqs. (5|6), by using the eight vertices as input. The value at the center of 



each face is assigned in the same way, but with the sum calculated over the four vertices 
corresponding to each face. Finally, the midpoint values of each of the twelve edges are 
calculated with the sum calculated over the vertices at the end points of the edges. The 
first iteration of this algorithm results in 27 subcubes. These steps, except the initial seeds 
of the eight vertices, are iteratively repeated for each of the 27 subcubes. Eventually, the 
number of subcubes will become (3-') d , where j is the iteration number and d = 3. Further 
details about this construction can be found in [33] . 

By using the detrending moving average (DMA) algorithm [51] . the Hurst exponent of 
the fractal structure can be subsequently estimated. The core of the DMA algorithm is the 
generalized variance o^ MA {s), that for d = 3 writes: 



DMA\ 



l,n2,n 3 ( r )\ , (7) 



where /ff(r) = fnih, h) is the fractional Brownian field with i\ = 1,2,..., N%, 12 = 
1, 2, N 2 and z 3 = 1, 2, A^ 3 . The function f ni ^ n2>n3 { r ) is given by 

fm,n2,n 3 {r) = ~ /ff(»"l - fci.ia - h,i3 ~ h), (8) 



v 

fci k 2 k s 



with the size of the subcubes (ni,n2,n 3 ) ranging from (3,3,3) to the maximum values 
( n imax, n 2max, n 3max) ■ v = ^1^-2^3 is the volume of the subcubes. The quantity V = (Ni — 
n imax)(N 2 — n 2ma x)(N 3 — n 3max ) is the volume of the fractal cube over the average / is 
defined. Eqs. ^ and fl8J) are defined for any geometry of the subarrays. In practice, it is 
computationally more suitable to use n\ = n 2 = n 3 to avoid spurious directionality and 
biases in the calculations. In Fig. [2j the log-log plots of (Tdma{ s ) vs s are shown f° r fractal 
cubes generated according to the above described procedure. The cubes have Hurst exponent 
H = 0.1, H = 0.2, H = 0.3, H = 0.4 and H = 0.5 respectively. The log-log plots of a 2 DMA (s) 
as a function of s are straight lines according to the power-law behavior: 

°Ima(s) « {n\ + n 2 2 + nl) H ex s H , (9) 

because of the fundamental property of fractional Brownian functions. 

By using these algorithms, the fractality of the snow structure can be related to the 
snow density by mapping the fractional Brownian field ///(r) to a density field p(r). In this 
framework, the Hurst exponent, varying as a continuous parameter, should be intended as an 
index of specific snow compactness. Different snow textures have been simulated by varying 
the minimum value of the density p m in between and 917kg/m 3 , while the maximum density 
is constant and equal to the ice density p max = p ice = 917 kg/m 3 . 

In Fig. [3] snow structures corresponding to cubes with size r$ = 100 mm and granular size 
t 'grain = 0.25 mm, with H = 0.1 are shown. The density ranges from 900 to 917 kg/ni (a) 
and from to 917 kg/m 3 (b). One can observe the difference between the more compact solid 
ice structure (a) and the almost fully porous media featured by several areas of lower density 
(b). Colors are scaled in such a way that darker areas correspond to higher densities. Finally, 
the average density peerage of the fractals generated according to the above procedure has 
been calculated. The results are plotted in Fig. |4j One can notice that the average density 
decreases more rapidly with lower values of p m in as the Hurst exponent increases, while 
the average density is practically unchanged as H is changed, by taking p min close to the 
value 917 kg/m 3 . The present approach might have interesting applications for monitoring 
ice losses and snow metamorphism. By independent measures of snow density and Hurst 
exponent, one can map the morphological evolution of snow/ice by using curves similar to 
those of Fig. |] 
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IV. CONCLUSIONS 



Physical and mechanical properties of snow are usually defined in terms of the specific 
density thanks to its simplicity of in-situ measurements. Unfortunately, the density is not 
univocally related to the snow microstructure, since different snow microstructure might 
exhibit the same global density. A scale invariant parameter is needed to quantify snow 
metamorphism in terms of the multiscale properties of snow density and porosity. We 
have proposed a fractal model for snow density based on (i) a generalized Menger sponge 
and (ii) a stochastic fractional Brownian field. The present approach shows that different 
Hurst exponents correspond to the same value of density, implying that density alone does 
not yield complete information about snow microstructure. Nonetheless, thanks to this 
model, one should be able to investigate how the local structure evolves according to the 
fractal dimension in relation to other physical properties. The present work is the first 
step towards the investigation of the scaling properties of snow in a fully three-dimensional 
fractal framework, relevant to the validation of experimental results, such as those reported 
in [TT| [TS]. and the description of the physical and mechanical properties which are of great 
interest for many application areas. 
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9 

(a) (b) (c) 

FIG. 1: Solid homogeneous cube at the initial stage of the iteration process (a), Menger sponge 
at the first iteration (b), Menger sponge at the second iteration (c). Size, density and porosity 
are respectively: r Q , p Q and <j) = (a); r\ = 3r c , pi = (20/27)p o and 4>\ = 7/27 (b); = 9r Q , 
p 2 = (400/729)/) o and <f> 2 = 329/729 (c) [H]. 
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FIG. 2: Log- log plot of the function a dm A as a function of scale s. The Hurst exponent is 
estimated by the slope of the best fit (red lines), respectively H = 0.1, H = 0.2, H = 0.3, H = 0.4 
and H = 0.5. 
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(a) (b) 

FIG. 3: Fractal cubes generated according to the three-dimensional fractional Brownian model 
presented in Section III. The Hurst exponent is H = 0.1 and the density p{r) ranges respectively 
between 800 — 917 kg/m 3 (a) and — 917 kg/m 3 (b). The ratio between the cube edge and the 
grain size is 400, implying that there are 400 x 400 x 400 values in each cube. 
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FIG. 4: Snow density as a function of the Hurst exponent. The different curves have been 
obtained by generating fractional Brownian functions defined over a cubic lattice [35J. Then, 
the fractional function has been mapped to a density function. In order to simulate different 
snow microstructures, the minimum value of the density p m in has been varied between and 917 
kg/m 3 , while the maximum value of the density is kept constant and equal to the ice density 
Pmax = Pice = 917 kg/m 3 . Different curves, from bottom to top, correspond to values of p m in 
ranging from to 917 kg/m 3 , with step 100 kg/m 3 . Dashed horizontal lines indicate density 
ranges for dry snow (50 < p < 200); snow (200 < p < 550); firn (550 < p < 820) and porous ice 
(820 < p < 917) after g|. 
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TABLE I: Classification of snow and ice (first column) according to the density (second column). 
The values of the Hurst exponent H have been calculated for a sample with linear size r snow = 
100 mm for different grain sizes r gra in on the basis of the generalized Menger sponge model discussed 
in Section II. The snow type and density are taken from the classification reported in [3]. 



Snow type 


Density [kg/m 3 ] 


r grain = 0.05 mm 


rgrain = 0.15 mm 


r grain = 0.25 mm 


Dry snow 


50 < p < 200 


0.3827 < H < 0.2003 


0.4474 < H < 0.2342 


0.4855 < H < 0.2542 


Snow 


200 < p < 550 


0.2003 < H < 0.0673 


0.2342 < H < 0.0786 


0.2542 < H < 0.0853 


Firn 


550 < p < 820 


0.0673 < H < 0.0147 


0.0786 < H < 0.0172 


0.0853 < H < 0.0187 


Porous ice 


820 < p < 917 


0.0147 < H < 0.0000 


0.0172 < H < 0.0000 


0.0187 < H < 0.0000 
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